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Introduction 


The Generalized Fluid System Simulation Package [2] (GFSSP) is a general-purpose computer 
program developed at NASA/MSFC for analyzing steady state and transient flow rates, 
pressures, temperatures, and concentrations in a complex flow network. The code can handle 
compressible and incompressible flows as well as phase change and mixture thermodynamics. 
The program which was developed out of need for an easy to use system level simulation tool for 
complex flow networks, has been used for the following purposes to name a few: Space Shuttle 
Main Engine (SSME) High Pressure Oxidizer Turbopump Secondary Flow Circuits, Axial 
Thrust Balance of the Fastrac Engine Turbopump, Pressurized Propellant Feed System for the 
Propulsion Test Article at Stennis Space Center, X-34 Main Propulsion System, X-33 Reaction 
Control System and Thermal Protection System, and International Space Station Environmental 
Control and Life Support System design. There has been an increasing demand for implementing 
a combustion simulation capability into GFSSP in order to increase its system level simulation 
capability of a propulsion system starting from the propellant tanks up to the thruster nozzle for 
spacecraft as well as launch vehicles. 

The present work was undertaken for addressing this need. In what follows, the model used for 
combustion of liquid hydrogen (LH 2 ) with liquid oxygen (LOX) using chemical equilibrium 
assumption, and the novel computational method developed for determining the equilibrium 
composition and temperature of the combustion products by application of the first and second 
laws of thermodynamics will be described. The modular FORTRAN code developed as a 
subroutine that can be incorporated into any flow network code with little effort has been 
successfully implemented in GFSSP as the preliminary runs indicate. The code provides 
capability of modeling the heat transfer rate to the coolants for parametric analysis in system 
design. 

Theoretical Model 


The current investigation considers the problem of combustion of LH 2 with LOX in a combustor 
with the assumption of H 2 0, II 2 , 0 2 , OH, H, and O constituting the products of 
combustion. The schematic diagram of the problem is depicted in Figure 1. 

x lb mol //., 


y lb mol 0 2 


Figure 1: Schematic diagram of the combustion chamber. 

The problem is solved by using chemical equilibrium assumption for steady flow case by 
applying the first and second laws of thermodynamics [4], An equivalent approach of using 
chemical equilibrium is the minimization of Gibbs free energy, which has been used in codes 
like CEA [1,5]. For given inlet conditions and exit pressure, the resulting equations, which turn 
out to be coupled, non-linear, algebraic equations, are solved simultaneously. Three cases are 
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considered; (a) stoichiometric case, (b) fuel rich case, and (c) oxygen rich case, for all possible 
environments. 


Stoichiometric Case : 

The chemistry mechanism used for this case is given by a set of reactions as follows 

2 H 2 + 0 2 o 2 H 2 0 
H 2 <=> 2 H 
0 2 +=> 20 

2 H 2 0 o H 2 + 2 OH 

and the changes in each specie in the form of fractions a , b , c , and d can be expressed as 

2 H 2 + 0 2 <=> 2 H 2 0 

-2a -a +2 a 

H 2 <=> 2 H 
-b +2 b 
0 2 <=> 2(9 

-c +2c 

2 H 2 0 <=> H 2 + 2 OH 
— 2d + d + 2t/ 


which result in the following overall reaction 

2 H 2 + 0 2 -> (2a-2d)H 2 0 + (2-2 ci-b + d)H 2 + (1 — a — c)0 2 

(2d)OH + (2b) H + (2c) O 


The individual mole fractions of species in the products can be expressed in terms of fractions a , 
b , c , and d in the following form 


_ 2a - 2d 
2 ° 3-a + b + c + d 


_ 2-2 a-b + d 
2 3-n+h + c + r/ 


1 - a -c 
3-a + b + c + d 


Z OH 


2d 

3-a+b+c+d 


2b 

3-a+b+c+d 


2c 

3-a+b+c+d 


and the application of second law with chemical equilibrium assumption yields 


z 2 

z // 2 0 

t p\ 

z 2 

K - H 

f P 1 

z 2 -z 
^h 2 ^o 2 

KPo) 

z h 2 
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*3=^ 

Z 0, 


f nA 


\ P oJ 


K 


Z H 2 ' Z OH 


( 


4 z 2 


H.O 


Oj 


and the first law of thermodynamics in the absence of work is expressed as 

Qc.v. + (^H 2 ) r (^// 2 )r + (^0 2 )r (^0 2 )/• = (^H 2 o)p f Hl0 + ^H 2 o) p + (^tf 2 ) p i^H 2 ) p 

+ («o 2 ) p (^o 2 )p + ) p (hf OH + ^h 0H ) p 

+ (« w )p (hf H +Ah H ) p +{h 0 ) p {h° fg +Ah 0 ) p 

where the number of moles per unit time, h i , for any specie i , can also be expressed in terms of 
the fractions a , b , c , and d can be expressed as 


p. , (2ci -2d) 

( n H 2 o)p= r (Or 


p. , (2-2 a-b + d) f . , 

(«// 2 )p = r K,), 


(1-a-c) 

K 2 )p = — - — KA 


(«o//)p =y-(«/w)r 


( j l H)p = 7 J«,)r 


(n°) p =yO*/«)r 


where (/i to , ),. = ),. + (< ),. , « ), 


m , 


MW U 


(Or 


m. 


MW n 


and 


, . . (3-a + b + c + J).. . 

(*«W ) p = ( n tot ) r • 

The five equations to be solved simultaneously can be written in the following fonn after some 
manipulations 


f = 

fi = 

f 3 = 
/ 4 = 


\ P oJ 


(2a-2d) 2 (2-a + b + c + d) f ^ 
(2-2 a-b + d) 2 (1-a-c) 

(2 bf 

(3 - a + b + c + d) (2 - 2a - b + d) 


-i 


K, =0 


f p \ 


\ P 0 J 


■k 2 = 0 


(2c) 2 


r p \ 


(3 - a + b + c + d) (1 - a - c) 
(2-2a - b + d)(2d) 


\ P o J 


2 ^ P ^ 


(2a - 2<f)“ (3 — a + b + c + J) y 


if 3 =0 


■ if 4 — 0 


( 1 ) 

( 2 ) 

( 3 ) 

( 4 ) 
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( 5 ) 


fs ~ Qc.v. ("ff 2 ) r (^n 2 ) r + ("o 2 )r (A^0 2 )r ("ff 2 o) p ^ f„ l0 p 

— ("ff 2 ) p (A^H 2 )p ~ («0 2 ) p (^0 2 )p — ("off ) p (^/ 0 „ ^h(3H ) p 

- (77; + A/ V ) ;; - (i* 0 ), (A; o + Ah 0 ) p = 0 

The corresponding equations have been developed for the fuel rich case and oxygen rich case. 


Fuel rich case (FRC) 

The chemistry mechanism for this case uses the following reactions 

x// 2 + y0 2 <=> 2 yH 2 0 + (x-2 y)H 2 for x > 2y 

H 2 o 2 H 
0 2 <=> 20 

2 H 2 0 ^ H 2 + 2 OH 

and the changes in each specie in the form of fractions a , b , c , and d can be expressed as 

xH 2 + y0 2 <=> 2 yH 2 0 + (x-2 y)H 2 

— 2ya -ya +2 ya +(x-2 y)a 

H 2 <=> 2 H 
-b +2 b 

0 2 <=> 20 
-c +2 c 

2 H 2 0 <=> H 2 + 2 OH 

— 2d + d + 2d 

The overall reaction is given by 

xH 2 + y0 2 — > ( 2ya-2d)H 2 0 + [x(l + a)-4ya -b + d]H 2 

+ (y- ya-c)0 2 + (2 d)OH + (2b) H + (2c) O 


Individual mole fractions are given as 

(2 ya - 2d) 


Z HjO ~ 


x(l + a) + y — 3 ya + b + c + d 

y -ya-c 

x(l + a) + y - 3 ya + b + c + d 

2b 


Z H = 


Z H 2 = 


-OH 


z o = 


x(l + a) - 4 ya - b + d 
x(l + a)+ v - 3 ya + b + c + d 

2d 

x(l + a) + y - 3 ya + b + c + d 

2c 


x(l + a) + y - 3 ya + b + c + d x(l + a) + y — 3ya + b + c + d 

With equilibrium constants given in terms of the mole fractions of products 
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*1 = 


7 2 >’ 

'//,0 


( n V V 


z 2y ■ z y 


\ P oJ 


K,=^L 


( rt \ 




\ p „j 


k,=Zl 


( rt \ 


y p oj 


k 4 = 


Z H , ' Z OH 


( T^\ 


■ H-,0 


\ P o ) 


The five coupled, non-linear algebraic equations for this case become 


f T,\~ y 


f 

fi 

h 

h 


K P oJ 


(2 ya - 2 d) 2y [x(l + a) + y - 3 ya + b + c + d]' 

[x(l + a) - 4 ya -b + d ] 2v (y - ya - c) y 

(2 Vf 

[x(l + a) + y - 3 ya + b + c + r/][x(l + a) - 4 ya — b + d\ 

. (2 cf 

[x(l + a) + y- 3 ya + b + c + d\y - ya - c ) 

[x(l + a) - 4 ya -b + d ](2 d) 


-K,= 0 


f P " 


r p\ 

k~P 0 j 


\ P o J 


-^3=0 


-^2=0 


r p \ 


(2ya-2d)° [x(l + a) + v -3ya + b + c + d]yP a j 


-K 4 = 0 


fs Qc.v. if l H 2 )/• (Ah Hl )/• («0 2 ) r (^h 0 ) r P (hf„ 20 + A h H0 ) p 

— (Mh 2 ) p (^b lr2 ) p — (h 02 ) p (Ah 02 ) p — (bon ) P (hf OH ^h 0H ) p 
-( n H ) P {h° fH +A h H ) p -{h 0 ) p {h° fo +Ah o ) p =0 


( 1 *) 

(2*) 

( 3 *) 

(4*) 

(5*) 


where the number of moles per unit time, h i , for any specie i, are expressed in terms of the 
fractions a , b , c , and d accordingly. 


Oxygen rich case (ORC) 

The chemistry mechanism for this case uses the following reactions 

xH 2 + y0 2 <=> xH 2 0 + (y~j)0 2 for x<2y 

H 2 o 2 H 

0 2 <=> 20 

2H 2 0 <=> H 2 + 2 OH 

and the changes in each specie in the form of fractions a , b , c , and d can be expressed as 
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xH 2 + y 0 2 

<=> 

xH 2 0 

+ (y-|)o 2 

X 

-xa a 

2 


+ xa 

+ (T- 


H 2 

<=> 2 H 



b 

+ 2 b 



0 2 

<=> 20 



-c 

+ 2c 


2 H 2 0 


H 2 + 

2 OH 

-2d 


+ d 

+ 2d 


The overall reaction is given by 

xH 2 + y0 2 — > (xa-2d)H 1 0 + [x(l — a) — b + d]H 2 

+ [y-(y — x)a-c\0 2 + {2d) OH + (2b) H + (2c) O 


Individual mole fractions are given as 

(xa - 2d) 


x(\ - a)-b + d 


■h 2 o 


z o , = 


x + v + (y - x)a + b + c + d 

y + (y - x)a - c 
x + y + (y - x)a + b + c + d 

2b 

x + y + (y - x)a + b + c + d 


Z OH ~ 


x + y + (y - x)a + b + c + d 

2 d 

x + y + (y - x)a + b + c + d 

2c 

x + y + (y - x)a + b + c + d 


With equilibrium constants given in terms of the mole fractions of products 


Kr = 


" h 2 o 


x/2 


' P 'l 

-x/2 . 

Z„ 

( p ' 


k 2 = h 


l P o) 

z h 2 

K P o) 




( 


K P oJ 


k 4 = 


Z H , ' Z OH 


( rt \ 




K P oJ 


The five coupled, non-linear algebraic equations for this case become 
(xa - 2 d) x [x + v + (y - x)a + b + c + d]' 


f 

f 2 


x/2 f p\~ xl2 


[x(l - a) - b + d]' [y + (y - x)a - c\ xn 

(2 bf 

[x + y + (y - x)a + b + c + d]x(l - a) —b + d\ 


K P oJ 


-K,= 0 


f p \ 


\ P o J 


-K 2 = o 


(1**) 

( 2 **) 
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A 

A 


(2c) 2 


f p \ 


[ x + v + (y - x)a +b + c + d~\ [y + (y- x)a - c ] 
[x(\-a)-b + d](2d) 2 


K P oJ 


-K 3 = 0 


r P \ 


(. xa - 2d) 2 [x + v + ( v - x)a + b + c + d] 




~K,=0 


(3**) 

(5**) 


f 5 Qc.v. (*G 2 )r i^H 2 )/• («0 2 )r (^0 2 )r i^ Hl o) p ^ f Hl o ^H 2 o) p 

— (^// 2 ) p (Ah Hl )p - («0 2 )p C^0 2 )p ~(Aoh) p (.hf OH + £dl OH ) p 

- («*), (1; + aa*), - (i* 0 ), (A; o + a/^ = 0 

where the number of moles per unit time, h l , for any specie i , are expressed in terms of the 
fractions a , b , c , and d accordingly. 


Solution Method 


It should be noted that whichever of the three cases is under consideration, the set of equations to 
be solved will be equations 1-5, or l*-5*, or 1 **_ 5 **. i n these equations the chemical 
equilibrium constants K x , K 2 , K , and A" 4 as well as all the A h i terms for the reactants and 
products are all functions of temperature. These functions have been correlated with least squares 
method with third or fourth order polynomials in three ranges of temperature from 0°R up to 
10000°R . 

Newton-Raphson Solution Scheme: 

The application of the Newton-Raphson method involves the following 7 steps [3]: 

1 . Develop the governing equations. The equations have to be in the following form: 


/l(x 1; x 2 , x 3 , x 4 , x 5 ) = 0 

(6a) 

X 

C/i 

II 

o 

(6b) 

/ 3 (x,, x 2 , x 3 , x 4 , x 5 ) = 0 

(6c) 

/ 4 (x,, x 2 , x 3 , x 4 , x 5 ) = 0 

(6d) 

/ 5 (x p x 2 , x 3 ,x 4 ,x 5 ) = 0 

(6e) 


Note that equations 1-5, l*-5*, and 1 **-5** are already in this form with a, b, c, d and 
T corresponding to x l5 x 2 , x 3 , x 4 , and x 5 respectively. 

2. Guess a solution for the equations: 

Guess x* , x*, x*, x*, x* as an initial solution for the governing equations. 

3. Calculate the residual of each equation. 

When the guessed solutions are substituted into equations (6a-6e), the right hand sides of the 
equations, which are not zero, are the residuals. 

/ (x* , x* , x* , x* , x* ) = R, (7a) 
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The intent of the solution scheme is to correct x* , x* , x* , x* , and x* with a set of 
?! ?! ?! ?! ?! 

corrections x t , x 2 , x 3 , x 4 , and x 5 such that R,, R 2 , R } , R 4 , and R 5 are zero. 

4. Develop a set of correction equations for all variables. 

First construct the matrix of influence coefficient: 


¥ 

¥ 

¥ 

5/ 

5/ 

dx l 

dx 2 

dx 3 

5x 4 

dx 5 

¥i 

¥2 

5/ 

5/ 

5/ 

dx l 

dx 2 

dx 3 

5x 4 

dx 5 

¥3 

¥3 

5/ 

5/ 

5/ 

dx, 

dx 2 

dx 3 

5x 4 

5x 5 

¥4 

¥4 

5/4 

5/4 

5/4 

dx. 

dx 2 

dx 3 

5x 4 

dx 5 

5/s 

5/s 

5/s 

5/s 

5/s 

cbq 

c)x 2 

fix'. 

5x 4 

dx 5 


Then construct the set of simultaneous equations for corrections: 


¥ 

5/ 

5/ 

5/ 

5/ 

dx. 

5x 2 

dx 3 

5x 4 

5x 5 

5/, 

5/, 

5/, 

5/z 

5/, 

dxj 

5x 2 

dx 3 

5x 4 

5x 5 

5/ 

5/ 

5/ 

5/3 

5/3 

dx t 

5x 2 

dx 3 

5x 4 

5x 5 

5/ 4 

5/4 

5/ 4 

5/i 

5/4 

3xj 

5x 2 

dx 3 

5x 4 

5x 5 

5/s 

5/s 

5/s 

5/s 

5/s 

dx 4 

5x 2 

dx 3 

5x 4 

5x 5 

e for 

?! 

Xj, J 

?! ?! 

c 2 , x 3 

?! 

, * 4 , 

and 


xj' 


¥ 

x 2 


r 2 

x 3 

> = < 

¥ 

x 4 


¥ 



¥, 


6. Apply correction to each variable. 

7. Iterate until corrections become smaller than a prescribed value. 
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Subroutines EQNS, SIMUL, PARDIF, and GAUSSY as described in reference [3] were 
modified for the problem under investigation. EQNS provides the simultaneous equations to be 
solved; SIMUL performs the Newton-Raphson solution of equations in EQNS with application 
of PARDIF and GAUSSY appropriately; PARDIF extracts the partial derivatives; and GAUSSY 
solves the resulting set of linear simultaneous equations by Gauss elimination method. 

Conclusion 


A modular FORTRAN code was developed for applying the first and second laws of 
thennodynamics to the combustion of LH 2 with LOX. Chemical equilibrium assumption based 
on the minimization of Gibbs free energy was employed. Case studies were performed with 
stoichiometric, fuel-rich and oxygen-rich cases. All these yielded very good results compared 
with hand calculations. Preliminary runs also indicate successful integration of this modular code 
with GFSSP. This modular code enhances the capability of GFSSP in performing system level 
simulation of a complete propulsion system using LH 2 /LOX, starting from the propellant tanks 
up to the thruster nozzle. 
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